Análisis de Datos I
Tarea 5

Librerías.

library(datasets)
library(FactoMineR) 
library(factoextra)
library(MASS)

Ejercicio 1.

Para la tabla de datos eurodist que viene con el paquete datasets el cual contiene las distancias entre algunas de las ciudades más importantes de Europa ejecute un Escalamiento Multidimensional y luego compare el resultado con el mapa de Europa.

Primero se obtienen los datos de las distancias apartir del paquete datasets.

matriz_distancias <- eurodist

Una vez que se cuenta con las distancias, se ejecuta el ACP o Multidimensional Scalling, a partir de la matriz de distancias.

MDS <- cmdscale(matriz_distancias,eig=TRUE, k=2)
MDS

## $points
##                         [,1]        [,2]
## Athens           2290.274680  1798.80293
## Barcelona        -825.382790   546.81148
## Brussels           59.183341  -367.08135
## Calais            -82.845973  -429.91466
## Cherbourg        -352.499435  -290.90843
## Cologne           293.689633  -405.31194
## Copenhagen        681.931545 -1108.64478
## Geneva             -9.423364   240.40600
## Gibraltar       -2048.449113   642.45854
## Hamburg           561.108970  -773.36929
## Hook of Holland   164.921799  -549.36704
## Lisbon          -1935.040811    49.12514
## Lyons            -226.423236   187.08779
## Madrid          -1423.353697   305.87513
## Marseilles       -299.498710   388.80726
## Milan             260.878046   416.67381
## Munich            587.675679    81.18224
## Paris            -156.836257  -211.13911
## Rome              709.413282  1109.36665
## Stockholm         839.445911 -1836.79055
## Vienna            911.230500   205.93020
## 
## $eig
##  [1]  1.953838e+07  1.185656e+07  1.528844e+06  1.118742e+06  7.893472e+05
##  [6]  5.816552e+05  2.623192e+05  1.925976e+05  1.450845e+05  1.079673e+05
## [11]  5.139484e+04 -3.259629e-09 -9.496124e+03 -5.305820e+04 -1.322166e+05
## [16] -2.573360e+05 -3.326719e+05 -5.162523e+05 -9.191491e+05 -1.006504e+06
## [21] -2.251844e+06
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 0.7537543 0.8679134

Finalmente se gráfica a los individuos y se muestra una comparación con el mapa verdadero.

Ubicaciones apartir del MDS:

x <- MDS$points[,1]
y <- -1*MDS$points[,2]
plot(x, y, xlab="Componente 1", ylab="Componente 2",main="MDS",pch = 19)
text(x, y, labels = labels(matriz_distancias), cex=0.85 , pos = 1)

Mapa real:

Como es posible observar, los resultados del MDS son muy similares al mapa real, con tan solo un para de diferencias que son atribuibles a no tener las distancias exactas o por lo contrario, al mapa con el que se cuenta.

Ejercicio 2.

Programe el algoritmo para Escalamiento Multidimensional (Mutidimensional Scaling) visto en clase y con la tabla de estudiantes compare con cmdscale() de R y veri que los resultados con un ejemplo.

A continuación, se presenta una función que realiza el MDS paso a paso.

MDS_algoritmo <- function(X, distancias) {
  
  # Paso 1: Con las distancias, se obtiene la matriz B 
  n <- nrow(X) # cantidad de individuos
  B_matriz <- matrix(NA, n, n)
  suma3 <- sum(distancias^2)
  
  for (r in 1: n) {
    suma1 <- 0
    suma2 <- 0
    for (s in 1: n ) {
      suma1 <- sum((distancias[,s])^2)
      suma2 <- sum((distancias[r,])^2)
      B_matriz[r,s] <- -(1/2)*(distancias[r,s]^2 - (1/n)*(suma1 + suma2) 
                               + (1/n^2)*suma3) 
    }
  }
  
  # Paso 2 : Se calculan los valores y vectores propios de B.
  B_e <- eigen(B_matriz)
  valores_propios <- B_e$values
  
  # Los valores propios muy pequeños cercanos a cero, se ignoran para el cálculo 
  # de las coordenadas. Se define un umbral para los valores propios 
  umbral <- 1e-10 
  
  # Se filtran los valores propios mayores que el umbral
  valores_propios_filtrados <- valores_propios[valores_propios > umbral]
  # Se obtienen los vectores propios asociados a los valores propios filtrados
  m <- length(valores_propios_filtrados)
  vectores_propios <- B_e$vectors[, 1:m]
  
  # Paso 3 : Se calculan las coordenadas
  coordenadas <- matrix(NA, n, m)
  rownames(coordenadas) <- rownames(X)
  
  for(i in 1: n) {
    for(j in 1: m){
      coordenadas[i,j] <- (sqrt(valores_propios_filtrados[j])*
                             vectores_propios[i,j])
    }
  }
  
  # Paso 4: Graficar 
  x <- coordenadas[,1]
  y <- coordenadas[,2]
  
  grafico <- plot(x, y, xlab="Componente 1", ylab="Componente 2",
                  main="MDS algoritmo",pch = 19)
  text(x, y, labels = row.names(X), cex = 0.85 , pos = 1)
  
  resultado <- list("points" = coordenadas, "eig" = valores_propios, 
                    "plot"= grafico)
  
  return(resultado)
}

Comparación con cmdscales()

Primeramente, se carga la base de datos de estudiantes.

estudiantes_datos <- read.table('EjemploEstudiantes.csv', 
                                header=TRUE, sep=';',dec=',',row.names=1)

Se obtiene las distancias de cada par de individuos.

distancias <- dist(estudiantes_datos)
distancias_matriz <- as.matrix(distancias)

Se verifican los resultados obtenidos por ambas funciones

MDS_estudiantes <- cmdscale(distancias, k = 5, eig=TRUE)
MDS_estudiantes

## $points
##               [,1]       [,2]        [,3]         [,4]          [,5]
## Lucia   0.76471745  1.5817637  1.11186219  0.039908252 -0.0007926494
## Pedro  -1.66887794 -1.3919656  0.09067929 -0.053941739  0.1128975404
## Ines   -1.57822841 -0.2994960  0.48752985  0.552652786 -0.1377925976
## Luis    2.60701317 -1.3202040 -0.46230941  0.784004734  0.0524431575
## Andres  1.43877557  1.3356687 -0.67985389 -0.189277341 -0.1117444886
## Ana    -2.34790534 -0.3880845 -0.12895699  0.007628838 -0.0253542859
## Carlos  0.89372557  1.5189012 -0.38893244 -0.124247926 -0.0093618642
## Jose   -2.64984571 -0.4254636 -0.46447580 -0.316066094 -0.0162142654
## Sonia   2.62959083 -2.1833951  0.40705140 -0.658738298 -0.0317927290
## Maria  -0.08896518  1.5722752  0.02740580 -0.041923214  0.1677121823
## 
## $eig
##  [1]  3.498309e+01  1.793415e+01  2.708155e+00  1.511504e+00  7.710194e-02
##  [6]  3.775779e-15  2.673082e-15  2.327544e-15 -9.116253e-16 -3.093847e-15
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 1 1

x <- MDS_estudiantes$points[,1]
y <- MDS_estudiantes$points[,2]

plot(x, y, xlab="Componente 1", ylab="Componente 2",
                main="MDS",pch = 19)
text(x, y, labels = row.names(estudiantes_datos), cex = 0.85 , pos = 1)

MDS_algoritmo_estudiantes <- MDS_algoritmo(estudiantes_datos, distancias_matriz)

Las coordenadas de los individuos son iguales. Por otro lado, los primeros 5 valores propios dados por ambos son los mismos. En cuanto a los restantes 4, estos difieren,pues, al ser números cercanos a cero por cuestiones de cálculo de la función cmdscale la precisión númerica cambia en comparación con eigen.Por último, el plano principal es el mismo.

Ejercicio 3.

Pruebe la ecuación (1.3) de la filminas vistas en clase, la cual establece una fórmula para el cálculo de brs.

Ejercicio 4.

Pruebe el Teorema que establece la relación entre los valores y vectores propios del algoritmo anterior en relación a los valores y vectores propios del ACP.

Ejercicio 5.

Del documento adjunto MDS no métrico.pdf leer la sección MDS no métrico y replicar en R los dos ejemplos que ah se presentan.

Ejemplo 1.

nombres_fila_columna1 <- c("M", "B", "V", "S", "SS", "LC")

m <- matrix(c(
  0,   627, 351, 550, 488, 603,
  627, 0,   361, 1043,565, 1113,
  351, 361, 0,   567, 564, 954,
  550, 1043,567, 0,   971, 950,
  488, 565, 564, 971, 0,   713,
  603, 1113,954, 950, 713, 0), 
  nrow = 6, byrow = TRUE)

m <- as.data.frame(m)
rownames(m) <- nombres_fila_columna1
colnames(m) <- nombres_fila_columna1
m <- as.matrix(m)

Calculamos Q

Q <- matrix(0, nrow = 6, ncol = 6)

sum1 <- c()
i <- 0
for (s in 1:6) {
  for (r in 1:6) {
    i <- i + (m[r,s])^2
    sum1[s] <- i
  }
  i <- 0 
}
  
sum2 <- c()
j <- 0
for (r in 1:6) {
  for (s in 1:6) {
    j <- j + (m[r,s])^2
    sum2[r] <- j
  }
  j <- 0 
}

sum3 <- 0
for (r in 1:6) {
  for (s in 1:6) {
    sum3 <- sum3 + (m[r,s])^2
  }
}

for (r in 1:6) {
  for (s in 1:6) {
    Q[r,s] <- -(1/2)*(((m[r,s])^2) - (1/6)*sum1[s] - (1/6)*sum2[r] + 
                        (1/6^2)*sum3)
  }
}

Q

##           [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,]  11759.44  -39079.22  -17954.39   38559.11  -31804.89   38519.94
## [2,] -39079.22  303211.11  124211.44 -208389.56   73380.44 -253334.22
## [3,] -17954.39  124211.44   75532.78   60951.28  -39894.22 -202846.89
## [4,]  38559.11 -208389.56   60951.28  367858.78 -206103.72  -52875.89
## [5,] -31804.89   73380.44  -39894.22 -206103.72  162774.78   41647.61
## [6,]  38519.94 -253334.22 -202846.89  -52875.89   41647.61  428889.44

Encontramos los valores y vectores propios de Q

propios <- eigen(Q)
valores_propios <- propios$values
valores_propios <- valores_propios[-c(5,6)]
vectores_propios <- propios$vectors
vectores_propios <- vectores_propios[,-c(5,6)]
valores_propios

## [1] 737916.932 591058.808  59468.730   1035.536

vectores_propios

##            [,1]        [,2]         [,3]        [,4]
## [1,]  0.0959729 -0.04429354  0.256942195  0.85659994
## [2,] -0.6270080  0.14004992  0.415521275 -0.15933263
## [3,] -0.2832133 -0.25835017  0.009409856 -0.31296904
## [4,]  0.2934011 -0.72163136 -0.220480281 -0.12846516
## [5,] -0.1240990  0.44165813 -0.781233037  0.08849908
## [6,]  0.6449462  0.44256702  0.319839993 -0.34433219

Calculamos las coordenadas

coordenadas <- matrix(0, nrow = 6, ncol = 4)
for (i in 1:6) {
  for (j in 1:4) {
    coordenadas[i,j] <- sqrt(valores_propios[j])*vectores_propios[i,j]
  }
}

coordenadas[,1:2]

##            [,1]       [,2]
## [1,]   82.44273  -34.05303
## [2,] -538.61294  107.67087
## [3,] -243.28619 -198.62051
## [4,]  252.03769 -554.79271
## [5,] -106.60360  339.54831
## [6,]  554.02232  340.24707

Graficamos

plot(coordenadas[,1], coordenadas[,2], col = "orange", pch = 16,
     xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(coordenadas[,1], coordenadas[,2], 
     labels = rownames(m), pos = 4, cex = 0.7, col = "black")

Ejemplo 2.

nombres_fila_columna2 <- c("Atlanta", "Chicago", "Denver", "Houston", 
                           "L. Angeles", "Miami", "N York", 
                           "S Francisco", "Seattle", "Washington")
aire.dist <- matrix(c(0, 587, 1212, 701, 1936, 604, 748, 2139, 218, 543,
                      587, 0, 920, 940, 1745, 1188, 713, 1858, 1737, 597,
                      1212, 920, 0, 879, 831, 1726, 1631, 949, 1021, 1494,
                      701, 940, 879, 0, 1374, 968, 1420, 1645, 1891, 1220,
                      1936, 1745, 831, 1374, 0, 2339, 2451, 347, 959, 2300,
                      604, 1188, 1726, 968, 2339, 0, 1092, 2594, 2734, 923,
                      748, 713, 1631, 1420, 2451, 1092, 0, 2571, 2408, 205,
                      2139, 1858, 949, 1645, 347, 2594, 2571, 0, 678, 2442,
                      218, 1737, 1021, 1891, 959, 2734, 2408, 678, 0, 2329,
                      543, 597, 1494, 1220, 2300, 923, 205, 2442, 2329, 0), 
                    nrow = 10, byrow = TRUE)
aire.dist <- as.data.frame(aire.dist)
rownames(aire.dist) <- nombres_fila_columna2
colnames(aire.dist) <- nombres_fila_columna2

Se efectua un analisis clasico MDS metrico

aire.mds <- cmdscale(as.matrix(aire.dist),k=9,eig=T)

## Warning in cmdscale(as.matrix(aire.dist), k = 9, eig = T): only 5 of the first
## 9 eigenvalues are > 0

aire.mds

## $points
##                   [,1]       [,2]       [,3]         [,4]          [,5]
## Atlanta      -434.7588  724.22221  440.92520   0.18579147 -0.0125796337
## Chicago      -412.6102   55.04016 -370.93031   4.39607623 12.6754984940
## Denver        468.1952 -180.65789 -213.57331  30.40856601 -9.5854646599
## Houston      -175.5816 -515.22265  362.83981   9.48713270 -4.8603541375
## L. Angeles   1206.6772 -465.63705   56.52609   1.34143943  6.8086243199
## Miami       -1161.6875 -477.98261  479.59934 -13.79783476  2.2781800452
## N York      -1115.5609  199.79247 -429.66594 -29.39692855 -7.1366804644
## S Francisco  1422.6887 -308.65595 -205.51788 -26.06309834 -1.9830609621
## Seattle      1221.5351  887.20174  170.44892  -0.06998691 -0.0000894326
## Washington  -1018.8972   81.89956 -290.65190  23.50884273  1.8159264311
## 
## $eig
##  [1]  9.213705e+06  2.199924e+06  1.082863e+06  3.322361e+03  3.858824e+02
##  [6] -5.234966e-09 -9.323115e+01 -2.168535e+03 -9.090644e+03 -1.722963e+06
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 0.8781612 1.0000000

Se Calculan los autovalores

aire.mds$eig

##  [1]  9.213705e+06  2.199924e+06  1.082863e+06  3.322361e+03  3.858824e+02
##  [6] -5.234966e-09 -9.323115e+01 -2.168535e+03 -9.090644e+03 -1.722963e+06

Se normalizan los dos primeros autovalores

sum(abs(aire.mds$eig[1:2]))/sum(abs(aire.mds$eig))

## [1] 0.8018277

sum(aire.mds$eig[1:2]^2)/sum(aire.mds$eig^2)

## [1] 0.9558842

Se muestran las coordenadas de las ciudades en las dos dimensiones

aire.mds$points[,1:2]

##                   [,1]       [,2]
## Atlanta      -434.7588  724.22221
## Chicago      -412.6102   55.04016
## Denver        468.1952 -180.65789
## Houston      -175.5816 -515.22265
## L. Angeles   1206.6772 -465.63705
## Miami       -1161.6875 -477.98261
## N York      -1115.5609  199.79247
## S Francisco  1422.6887 -308.65595
## Seattle      1221.5351  887.20174
## Washington  -1018.8972   81.89956

Se dibujan las coordenadas de las ciudades en las dos dimensiones

par(pty="s")
plot(-aire.mds$points[,1],aire.mds$points[,2],type="n",xlab="Coordenada 1",
     ylab="Coordenada 2",xlim = c(-2000,1500),ylim=c(-2000,1500))
text(-aire.mds$points[,1],aire.mds$points[,2],labels=row.names(aire.dist))

otra manera de dibujar las coordenadas

plot(-aire.mds$points[,1], aire.mds$points[,2], col = "orange", pch = 16, 
     xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(-aire.mds$points[,1], aire.mds$points[,2], 
     labels = rownames(aire.mds$points), pos = 3, cex = 0.6, col = "black")

data(swiss)

El fichero original de datos viene recogido ya en MASS:

swiss

##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
## Broye             83.8        70.2          16         7    92.85
## Glane             92.4        67.8          14         8    97.16
## Gruyere           82.4        53.3          12         7    97.67
## Sarine            82.9        45.2          16        13    91.38
## Veveyse           87.1        64.5          14         6    98.61
## Aigle             64.1        62.0          21        12     8.52
## Aubonne           66.9        67.5          14         7     2.27
## Avenches          68.9        60.7          19        12     4.43
## Cossonay          61.7        69.3          22         5     2.82
## Echallens         68.3        72.6          18         2    24.20
## Grandson          71.7        34.0          17         8     3.30
## Lausanne          55.7        19.4          26        28    12.11
## La Vallee         54.3        15.2          31        20     2.15
## Lavaux            65.1        73.0          19         9     2.84
## Morges            65.5        59.8          22        10     5.23
## Moudon            65.0        55.1          14         3     4.52
## Nyone             56.6        50.9          22        12    15.14
## Orbe              57.4        54.1          20         6     4.20
## Oron              72.5        71.2          12         1     2.40
## Payerne           74.2        58.1          14         8     5.23
## Paysd'enhaut      72.0        63.5           6         3     2.56
## Rolle             60.5        60.8          16        10     7.72
## Vevey             58.3        26.8          25        19    18.46
## Yverdon           65.4        49.5          15         8     6.10
## Conthey           75.5        85.9           3         2    99.71
## Entremont         69.3        84.9           7         6    99.68
## Herens            77.3        89.7           5         2   100.00
## Martigwy          70.5        78.2          12         6    98.96
## Monthey           79.4        64.9           7         3    98.22
## St Maurice        65.0        75.9           9         9    99.06
## Sierre            92.2        84.6           3         3    99.46
## Sion              79.3        63.1          13        13    96.83
## Boudry            70.4        38.4          26        12     5.62
## La Chauxdfnd      65.7         7.7          29        11    13.79
## Le Locle          72.7        16.7          22        13    11.22
## Neuchatel         64.4        17.6          35        32    16.92
## Val de Ruz        77.6        37.6          15         7     4.97
## ValdeTravers      67.6        18.7          25         7     8.65
## V. De Geneve      35.0         1.2          37        53    42.34
## Rive Droite       44.7        46.6          16        29    50.43
## Rive Gauche       42.8        27.7          22        29    58.33
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6
## Broye                    23.6
## Glane                    24.9
## Gruyere                  21.0
## Sarine                   24.4
## Veveyse                  24.5
## Aigle                    16.5
## Aubonne                  19.1
## Avenches                 22.7
## Cossonay                 18.7
## Echallens                21.2
## Grandson                 20.0
## Lausanne                 20.2
## La Vallee                10.8
## Lavaux                   20.0
## Morges                   18.0
## Moudon                   22.4
## Nyone                    16.7
## Orbe                     15.3
## Oron                     21.0
## Payerne                  23.8
## Paysd'enhaut             18.0
## Rolle                    16.3
## Vevey                    20.9
## Yverdon                  22.5
## Conthey                  15.1
## Entremont                19.8
## Herens                   18.3
## Martigwy                 19.4
## Monthey                  20.2
## St Maurice               17.8
## Sierre                   16.3
## Sion                     18.1
## Boudry                   20.3
## La Chauxdfnd             20.5
## Le Locle                 18.9
## Neuchatel                23.0
## Val de Ruz               20.0
## ValdeTravers             19.5
## V. De Geneve             18.0
## Rive Droite              18.2
## Rive Gauche              19.3

swiss.x <- as.matrix(swiss[, -1])
swiss.dist <- dist(swiss.x)
swiss.mds <- isoMDS(swiss.dist)

## initial  value 2.979731 
## iter   5 value 2.431486
## iter  10 value 2.343353
## final  value 2.338839 
## converged

Las coordenadas son

swiss.mds$points

##                    [,1]        [,2]
## Courtelary    39.976648 -18.6652658
## Delemont     -42.144880 -15.8354737
## Franches-Mnt -49.183472 -23.5017064
## Moutier       10.640047  -7.7948873
## Neuveville    35.713299   4.2241084
## Porrentruy   -45.069549 -26.7640938
## Broye        -55.273717   3.1907784
## Glane        -58.641853   0.1769567
## Gruyere      -55.274386 -11.8305230
## Sarine       -46.800689 -18.2672079
## Veveyse      -59.271433  -2.7666474
## Aigle         27.639864  18.9592510
## Aubonne       31.354629  26.3662771
## Avenches      32.353450  19.3520700
## Cossonay      31.342613  28.5036282
## Echallens     10.799996  26.1967959
## Grandson      40.595874  -1.7999750
## Lausanne      36.972991 -25.1351295
## La Vallee     50.122185 -23.5248585
## Lavaux        30.346198  30.8619325
## Morges        31.803525  18.6695352
## Moudon        33.621794  16.2266437
## Nyone         25.423297   7.4959696
## Orbe          34.176435  14.5941476
## Oron          30.076709  32.2309096
## Payerne       31.559866  17.6683208
## Paysd'enhaut  32.631374  27.0802839
## Rolle         28.769137  19.1245372
## Vevey         29.599603 -16.5009238
## Yverdon       33.286505  10.5935022
## Conthey      -67.755640  18.1702425
## Entremont    -66.402167  14.4997009
## Herens       -68.766153  20.1878849
## Martigwy     -63.370151   8.9517543
## Monthey      -60.005854  -0.5353731
## St Maurice   -62.872702   7.3903698
## Sierre       -66.745414  16.5445662
## Sion         -56.775840  -4.1215356
## Boudry        37.427364  -1.3796084
## La Chauxdfnd  41.684937 -30.7139972
## Le Locle      39.469769 -20.2862591
## Neuchatel     34.177795 -33.7884464
## Val de Ruz    37.403309   1.2494697
## ValdeTravers  42.204515 -17.0228055
## V. De Geneve  17.524061 -65.0238668
## Rive Droite   -6.829233 -11.9126783
## Rive Gauche   -7.514655 -31.3383738

Se dibujan los puntos

plot(swiss.mds$points, type="n", xlab="Coordenada 1", ylab="Coordenada 2")
text(swiss.mds$points, labels=as.character(1:nrow(swiss.x)))

Otra manera de dibujar los puntos con los nombres de etiqueta

coordenadas <- swiss.mds$points
nombres <- rownames(swiss.mds$points)

plot(coordenadas[,1], coordenadas[,2], col = "orange", pch = 16, 
     xlab = "Coordenada X", ylab = "Coordenada Y", main = "Plot de Coordenadas")
text(coordenadas[,1], coordenadas[,2], labels = nombres, pos = 4, cex = 0.3, 
     col = "black")

Análisis de Datos I Tarea 5